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Abstract 

Analytic formulas are presented for three flavor neutrino oscillations in matter in the 
plane wave approximation. We calculate in particular the time evolution operator 
in both mass and flavor bases. We also find the transition probabilities expressed as 
functions of the vacuum mass squared differences, the vacuum mixing angles, and 
the matter density parameter. The application of this to neutrino oscillations for 
both atmospheric and long baseline neutrinos in a mantle-core-mantle step function 
model of the Earth's matter density profile is discussed. 

PACS: 14.60.Pq, 14.60.Lm, 13.15.+g, 96.40.Tv 



1 Introduction 



As an accumulating amount of data on neutrino oscillations is becoming ac- 
cessible, it is of interest to study three flavor neutrino oscillations. Here we 
would like to give analytic expressions for the neutrino oscillation probabilities 
in presence of matter expressed in the mixing matrix elements and the neu- 
trino energies or masses, i.e., incorporating the so called Mikheyev-Smirnov- 
Wolfenstein (MSW) effect [1,2]. These probability formulas are of interest for 
the solar neutrino problem, the atmospheric neutrino problem, and the long 
baseline (LBL) neutrino experiments. We here apply the formalism to neutri- 
nos traversing the Earth in a mantle-core-mantle step function model of the 
Earth's matter density profile. 
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We will assume that CP nonconservation is negligible at the present level of 
experimental accuracy [3]. Thus, the mixing matrix for the neutrinos is real. 

Previous work on models for three flavor neutrino oscillations in matter in- 
cludes works of Barger et al. [4], Kim and Sze [5], and Zaglauer and Schwarzer 
[6]. Our method is different from all these approaches and their parameter- 
izations also differ slightly from ours. In particular, we calculate the time 
evolution operator and do not use the auxiliary matter mixing angles. 

Approximate solutions for three flavor neutrino oscillations in matter have 
been presented by Kuo and Pantaleone [7] and Joshipura and Murthy [8]. Ap- 
proximate treatments have also been done by Toshev and Petcov [9]. D'Olivo 
and Oteo have made contributions by using an approximative Magnus expan- 
sion for the time evolution operator [10]. Extensive numerical investigations for 
matter enhanced three flavor oscillations have been made by Fogli et al. [11]. 
Studies of neutrino oscillations in the Earth has been performed by several 
authors [12-18]. 



2 The evolution operator 



Let the flavor state basis and mass eigenstate basis be given bjHf = { | v a ) } a =e„ 
and 7i rn = {|^ a )}a=i, respectively. Then the flavor states \v a ) G Tif can be 
obtained as superpositions of the mass eigenstates \v a ) G 7i m , or vice versa. 
The bases Hf and H m are of course just two different representations of the 
same Hilbert space H. 

In the present analysis, we will use the plane wave approximation to describe 
neutrino oscillations. In this approximation, a neutrino flavor state \v a ) is a 
linear combination of neutrino mass eigenstates |z/ a ) such that [19] 

k> = ££Ck>, (i) 

a=l 

where a = e, fi, r. In what follows, we will use the short-hand notations \a) = 
\u a ) and | a) = \u a ) for the flavor states and the mass eigenstates, respectively. 

For the components of a state in the flavor state basis and mass eigenstate 
basis, respectively, they are related to each other by 

i>s = Uij m , (2) 
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where 



i>S = (V>a) 



eH f and = (ip a ) 



^2 



A convenient parameterization for U = U(9i, 9 2 , 9 3 ) is given by [20] 



U = 



— S3C1 — S1S2C3 C\C% — S1S2S3 S±C 2 
S\S 3 — S2C1C3 —S1C3 — S2S3C1 C1C2 J 



(3) 



where Si = sin^ and C, = cos#j for i = 1,2,3. This is the standard repre- 
sentation of the mixing matrix. The quantities 9 iy where i = 1,2,3, are the 
vacuum mixing angles. Since we have put the CP phase equal to zero in the 
mixing matrix, this means that U* a = U aa for a = e, /i, r and a = 1, 2, 3. 

In the mass eigenstate basis, the Hamiltonian J$? for the propagation of the 
neutrinos in vacuum is diagonal and given by 



( E x X 

E 2 
^ E 3 J 



(4) 



where E a = \Jm 2 a + p 2 , a = 1, 2, 3, are the energies of the neutrino mass 
eigenstates \a), a — 1, 2, 3 with masses m a , a = 1, 2, 3. We will assume p to be 
the same for all mass eigenstates. 

When neutrinos propagate in matter, there is an additional term coming from 
the presence of electrons in matter [2]. This term is diagonal in the flavor state 
basis and is given by 



( A 0^ 



Vf 








(5) 



where 



A = ±V2G F N e ~ ±-^=G F — p 
V2 m N 



is the matter density parameter. Here G F is the Fermi weak coupling constant, 
iV e is the electron density, is the nucleon mass, and p is the matter density. 
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The sign depends on whether we deal with neutrinos (+) or antineutrinos (— ). 
We will assume that the electron density N e (or the matter density p) is 
constant throughout the matter in which the neutrinos are propagating. In 
the mass basis, this piece of the Hamiltonian is V m = U^VfU, where again 
U is the mixing matrix. 

The unitary transformation that leads from the initial state tpf(0) in flavor 
basis at time t — of production of the neutrino, to the state of the same 
neutrino ipf(t) at the detector at time t is given by the operator Uf(t) = 
Uf(t,0), where Uf(t 2 ,ti) is the time evolution operator from time t± to time 
t 2 in flavor basis. This operator can be formally written as Uf(t) = e~ tHft . 
When the neutrinos are propagating through vacuum, the Hamiltonian in 
flavor basis is Hf = UH m U~~ l . In this case, the exponentiation of Hf can 
be performed easily: Uf(t) = e~ ft = U e~ %Uvnt XJ~ x , and the result can be 
expressed in closed form. In the case when the neutrinos propagate through 
matter, the Hamiltonian is not diagonal in either the mass eigenstate basis or 
the flavor state basis, and we have to calculate the evolution operator Uf(t) 
or Uf(L) = e-^ L = Ue~ ijeinL U- 1 if we set t — L. 

To do so it is convenient to introduce the traceless matrix T defined by T ee 
M'm — (tr J^)//3. The trace of the Hamiltonian in the mass basis ffl m ee 
H m + U~ l VfU is tr J^ m = Ei + E 2 + E 3 + A, and the matrix T can then be 
written as 



T = (T ab ) 

' AU 2 el -\A+\{E 12 + E lz ) AU el U e2 AU el U e3 N 

AU el U e2 AU 2 e2 -\A + \(E 21 +E 2 z) AU e2 U e3 

K AU el U e3 AU e2 U e3 AU%-tA+l(E 31 + E 32 ) / 

(6) 

where E ab = E a — E b . Of the six quantities E ab , where a,b = 1,2,3 and 
a ^ b, only two are linearly independent, since the E ab s fulfill the relations 
E ba = —E ab and E\ 2 + E 23 + E 3 \ = 0.0 Using Eq. (6), the evolution operator 
in the mass basis can be written as 

U m {L) ee e~^ L = <Pe~ iLT , (7) 

where ee e -* itr -^V3_ 

The flavor states can be expressed either as linear combinations of the vacuum 

3 Later, we will use the usual (vacuum) mass squared differences Am^ and Am^ 
instead of E 2 i and E 32 , which are related to each other by Am 21 = 2E U E 2 \ and 
Am| 2 = 2E V E 32 , where E v is the neutrino energy. 
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mass eigenstates (A = 0) in the basis 7i m as in Eq. (2) or the matter mass 
eigenstates (A ^ 0) in the basis Hm- 111 the latter case, the corresponding 
components are related to each other as 

i>i = U M ij M , (8) 

where U M = U M (9f r ,9^ , 0f*) is the unitary mixing matrix for matter and 
9^, i = 1,2,3, are the (auxiliary) matter mixing angles. 

Combining these expressions for the flavor state components, one obtains the 
following relation between the two different sets of mass eigenstate components 

1p M = Rlpm, (9) 

where 

R={U M )' l U. 

The matrix R is, of course, a unitary matrix (even orthogonal, since U and 
U M are real). This means that the matter mixing matrix can be expressed in 
the vacuum mixing matrix as 

U M {9\\ 9™, 9™) = U(9 ± , 9 2 , 9 3 )R~\ (10) 

The relations between the different bases are depicted in the following diagram: 

U 

i) m e H m - Hj 3 ip f 




4>m e T~Cm 

From this diagram one readily obtains 

M'm = RJ^mR' 1 or Jf M = RTR- 1 + i(trJ^)/, 

where J#m is the Hamiltonian in matter, which is diagonal in the basis TIm- 

Due to the invariance of the trace, we have 

T M = J^ M - §(tr JT M )/ = RTR' 1 . (11) 

and T M is a diagonal matrix with elements A a , a = 1,2,3, the eigenvalues of 
T. This implies that 

e -iLT M = Re ~iLT R -l ^ 
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Now, Cayley-Hamilton's theorem implies that, since T is a 3 x 3 matrix, the 
infinite series defining e~ lLT can be written as a second order polynomial in T 
[21]: 



-iLT 



= oqI - iLTax - L 2 T 2 a 2 , 



(13) 



where ao, Q>i, and a 2 are coefficients to be determined. Inserting Eq. (13) into 
Eq. (12) and using Eq. (11) gives a linear system of three equations that will 
determine the coefficients a , a ly and a 2 : 



-iLXi 



a — iLXidi — L \{a 2 



e -iL\ 2 = ao - iL\ 2 ax - L 2 }? 2 a 2 , 



(14) 



-1LX3 



= ao — iL\ 3 ai — L \\a 2 



where A a , a = 1,2,3, are the diagonal elements of Tm, or, equivalently, the 
eigenvalues of T, i.e., the solutions to the characteristic equation 



A 3 + c 2 A 2 + ciA + c = 



with 



(15) 



c 2 = -trT = 0, 
ci = detTtrT -1 , 
Co = — det T. 



When the system of equations (14) is solved for the o^'s, we obtain 



= E 



i 



a=! 3A2 + Cl 

From Eqs. (7) and (19) we then obtain 



(\ 2 a + Cl )I + \ a T + T 2 



(16) 
(17) 
(18) 



The coefficients Co, ci, and c 2 are all real and the eigenvalues A a , a = 1,2,3, 
can be expressed in closed form in terms of these [21]. 



(19) 



ML) = e-^ L = 0e-* iT = 0]>> 



iLX a 



a=l 



3X1 + Cl 



(\ 2 a + Cl )I + \ a T + T 2 



(20) 



The evolution operator for the neutrinos in flavor basis is thus given by 



Uf(L) = e~ iJf f L = Ue-^^lf- 1 = 0^ e 



iLX a 



a=l 



3X1 + ci 



(\l + Cl )I + \ a f + f 2 \ 
(21) 
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where T = XJTXJ~ X . Formula (21) is the final expression for the evolution 
operator. It expresses the time (or L) evolution directly in term of the mass 
squared differences and the vacuum mixing angles without introducing the 
auxiliary matter mixing angles. 



3 Transition probabilities 

The probability amplitude is defined as 

A a/3 = ((3\U f (L)\a), a,P = e,n,T. (22) 
Inserting Eq. (21) into Eq. (22) gives 

A _ A V" „-iL\a + C l)<W + KTgp + (T 2 ) a/3 

^af3-<P2^ e o\2 7 " ' y Z6 > 

o=l 6A a C l 

where 5 a/3 is Kronecker's delta. Note that T a/3 = Tp a and (T 2 ) a/3 = (T 2 ) /3q . 
Inserting L = into Eq. (21) yields 

r y- (^q + gQjag + \ a T a/3 + (T 2 ) a(3 

= £ 3A^ • (M) 

Hence, the transition probabilities in matter are 

p i, |2 r A y' y' (^a + c l)(W + AqT a/3 + 

Pa/3 = |A*/3| =^-42,L ox 2 T ~ 

q=l 6=1 dA a + C l 

q<fe 

„ ( A 6 + c l)<W + A fe T a/3 + (T 2 ) a/3 . 2 _ 

x —r- sm x ab , a,(3 = e,n,r, 

3\ z b + ci 

(25) 

where x ab = (A a - \ b )L/2. 

From unitarity there are only three independent transition probabilities, since 
the other three can be obtained from them, i.e., from the equations 

P ee + P efl + P eT = 1, (26) 

P^e + P^ + P„r = 1, (27) 

P re + P rM + P rr = 1. (28) 

Note that P eM = P Me , P er = P re , and P^ T = P Tii . If we choose P eM , P er , and 
P UT as the three independent ones, we thus have for a ^ f3 
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^ = -4 EE 



a=l 6=1 

a<b 



3 3 




sin x ab . 



(29) 



4 Applications and discussion 

The main results of our analysis are given by the time evolution operator for 
the neutrinos when passing through matter with constant matter density in 
Eq. (21) and the expression for the transition probabilities in Eqs. (25) and 
(29), expressed as finite sums of elementary functions in the elements of T. 

In our treatment the auxiliary mixing angles in matter, Of 1 , play no indepen- 
dent role and are not really needed. 

When the neutrinos travel through a series of matter densities with matter 
density parameters Ai, . . . ,A n and thicknesses L ± , . . . , L n , the total evolution 
operator is simply given by 



where L = X]" = i A and Uf(Lj) is calculated for A = Ai. Equation (30) gives 
a simple and fast algorithm to obtain the total evolution operator instead of 
using numerical integration. 

As an application to show the usefulness of our derived formulas, we have 
calculated the transition probabilities P a p for neutrino oscillations in a mantle- 
core-mantle step function model simulating the Earth's matter density profile. 
Let R ~ 6371 km be the radius of the Earth and r ~ 3486 km be the radius 
of the core. The thickness of the mantle is then R — r~ 2885 km, with matter 
density parameter A\ ~ 1.70 • 10~ 13 eV (pi ~ 4.5 g/cm 3 ), whereas the matter 
density parameter of the core is A 2 — 4.35 • 10 -13 eV (p 2 — 11.5 g/cm 3 ). 

Neutrinos traversing the Earth towards a detector close to the surface of the 
Earth, pass through the matter densities Ai,A 2 ,A\ of thicknesses L\ ) L 2) Li 
where the distances L iy % = 1,2, are functions of the nadir angle h, where 
h = 180° — Z ; 9 Z being the zenith angle. As h varies from to 90°, the cord 
L = L(h) of the neutrino passage through the Earth becomes shorter and 
shorter. At an angle larger than ho = arcsin(r/i?) ~ 33.17°, the distance 
L 2 = 0, and the neutrinos no longer traverse the core. 

For < h < h the distances L\ and L 2 are given by 



n 



U f (L) = l[U f (L i ) = U f (L n )...U f (L 1 ), 



(30) 
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r \ 



2 



L X =R \ cosh-\j(-J -sin 2 /i), (31) 



L2=2jR V© _sin2/l ' (32) 

L = 2L X + L 2 . (33) 



For h <h< 90°: 



L = 2Rcosh. (34) 



The mass squared differences (AM 2 = Am 2 2 and Am 2 = Am^) and the 
vacuum mixing angles (9\, 9 2 , #3) used here are chosen to correspond to those 
obtained from analyses of various neutrino oscillation data. For our illustration 
we have taken the following parameter values 

AM 2 = 3.2 • l(T 3 eV 2 , Am 2 = 0, AM 2 /10, 9 X = 45°, 9 2 = 5°, 9 3 = 45°. 

The values of AM 2 and 9\ are governed by atmospheric neutrino data [22] and 
the values of Am 2 and #3 by solar neutrino data [23] . The value of 9 2 is below 
the CHOOZ upper bound, which is sin 2 29 2 = 0.10 [24]. These choices are the 
most optimistic ones for obtaining any effects in LBL experiments from the 
sub- leading Am 2 scale [25]. We should mention though, that these data are 
taken from two flavor model analyses. 

The probabilities corresponding to various situations in this scenario are il- 
lustrated in Figs. 1-4. These figures correspond to the physically measurable 
quantities. From the figures one can see that there are several resonance phe- 
nomena superimposed on each other. 

Figures 1 a) - 1 d) show the results for the case of h — 0, one mass squared 
difference equal to 3.2 • 10~ 3 eV 2 [22], and the other one equal to 0. In Fig. 1 
a) we have included for comparison the corresponding result of a two flavor 
model. Figures 1 b) - 1 d) show the probabilities P ee , P ejti , and P^ T , respectively, 
in the three flavor model. The other probabilities can be obtained by unitarity 
from these three using Eqs. (26) - (28). We observe that although the survival 
probability P ee is the same for both cases, the transition probability P efJi is only 
half of that in the two flavor model. In a two flavor model, only one transition 
probability is needed, since if P ee is given and P^ e = P e/Ji , the other two follow 
from P efX = 1 - P ee and P m = 1 - P e(1 = P ee . 

In Figs. 2 a) - d) the second mass squared difference is chosen to be 1/10 of 
the first one. This affects the disappearance rate for low energies in the three 
flavor case, as can be seen from Fig. 2 b). It also modifies the appearance rates 
P eti and P^ T in Figs. 2 c) and 2 d), respectively. Figure 2 a) is again the same 
as Fig. 1 a). 
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When the small mass squared difference is much smaller than the large mass 
squared difference, then the two flavor model coincides with the three flavor 
model. This is not the case when the small mass squared difference is of 
comparable size to the large mass squared difference. See Figs. 1 a), 1 b), 2 
a), and 2 b) for a comparison. 

Figure 3 shows a contour plot of P ee as a function of the matter density 
parameter of the core A COTC (as it changes from the value of the mantle A\ up 
to its full value A 2 ) and the neutrino energy E u . The dark areas correspond 
to regions of large conversion probability for the assumed parameter values. 
When the small mass squared difference decreases, the structure below the 
main conversion valley goes away. Furthermore, the conversion regions above 
the main conversion valley are obviously due to the core (either purely to the 
core or to interference effects between the core and the mantle). 

Similarly, Fig. 4 shows a contour plot of P ee as a function of the nadir angle 
h and the neutrino energy E v . As in the previous figure, the dark areas cor- 
respond to regions of large conversion probability for the assumed parameter 
values. At h — ho, the neutrinos no longer traverse the core, but only the 
mantle. Again, when the small mass squared difference drops, the oscillations 
in the bottom of the plot (sub-GeV region) go away and no conversion takes 
place, i.e., P ee ~ 1 in this region. In the limit when the nadir angle goes to 
90°, the survival probability P ee becomes 1, because the traveling path length 
L of the neutrinos approaches 0. 

The results and the applications presented here are relevant to both the atmo- 
spheric neutrino experiments (e.g. Super-Kamiokande) and the LBL neutrino 
experiments. The present treatment assumes that in the results of the atmo- 
spheric and LBL experiments, the scattering and absorption of the neutrinos 
in matter is taken care of in the Monte Carlo programs that are used to an- 
alyze the data from these experiments. In a realistic analysis, we should also 
introduce damping due to uncertainties in the energy E v of the neutrinos and 
their path length L. 

Examining some of the future LBL neutrino experiments, we find their nadir 
angles to be h ~ 86.7°, h ~ 86.7°, and h ~ 88.9° for CERN-NGS (L = 743 
km), MINOS (L = 732 km), and K2K (L = 252 km), respectively [25,26]. See 
Fig. 4 for how the survival probability P ee depends on the neutrino energy for 
these fixed nadir angles. Note that for these three experiments the neutrinos 
only traverse the mantle and not the core. 

There are, at the moment, no planned LBL experiments in which the neutrinos 
also traverse the core. The baseline has in that case to be longer than L ~ 
10670 km. 
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a) b) 




£ v ( eV ) £ v ( eV ) 



Fig. 1. Transition probabilities as functions of the neutrino energy E v for the param- 
eter values h = 0, 6>i = 63 = 45° (bimaximal mixing), 9 2 = 5°, AM 2 = 3.2-1CT 3 eV 2 , 
and Am 2 = 0. a) P ee (two flavors; parameters values: 6 = 5° and Am 2 = 3.2 • 10~ 3 
eV 2 ), b) P ee (three flavors), c) P e)l (three flavors), and d) P^ T (three flavors). 



a) b) 




£ v (eV) £ v (eV) 

Fig. 2. Transition probabilities as functions of the neutrino energy E u , using the 
same parameter values as for Fig. 1 except that Am 2 = AM 2 /W. a) P ee (two 
flavors; parameters values: 6 = 5° and Am 2 = 3.2- 10~ 3 eV 2 ), b) P ee (three flavors), 
c) P eil (three flavors), and d) P^ T (three flavors). 
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Fig. 3. The survival probability P ee as a function of the core matter density param- 
eter A core e [A\,A2\ and the neutrino energy E u , using the same parameter values 
as for Fig. 2. 
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